Transcriptome Profiling and Chlorophyll Metabolic Pathway Analysis Reveal the Response of Nitraria tangutorum to Increased Nitrogen

To identify genes that respond to increased nitrogen and assess the involvement of the chlorophyll metabolic pathway and associated regulatory mechanisms in these responses, Nitraria tangutorum seedlings were subjected to four nitrogen concentrations (N0, N6, N36, and N60: 0, 6, 36, and 60 mmol·L−1 nitrogen, respectively). The N. tangutorum seedling leaf transcriptome was analyzed by high-throughput sequencing (Illumina HiSeq 4000), and 332,420 transcripts and 276,423 unigenes were identified. The numbers of differentially expressed genes (DEGs) were 4052 in N0 vs. N6, 6181 in N0 vs. N36, and 3937 in N0 vs. N60. Comparing N0 and N6, N0 and N36, and N0 and N60, we found 1101, 2222, and 1234 annotated DEGs in 113, 121, and 114 metabolic pathways, respectively, classified in the Kyoto Encyclopedia of Genes and Genomes database. Metabolic pathways with considerable accumulation were involved mainly in anthocyanin biosynthesis, carotenoid biosynthesis, porphyrin and chlorophyll metabolism, flavonoid biosynthesis, and amino acid metabolism. N36 increased δ-amino levulinic acid synthesis and upregulated expression of the magnesium chelatase H subunit, which promoted chlorophyll a synthesis. Hence, N36 stimulated chlorophyll synthesis rather than heme synthesis. These findings enrich our understanding of the N. tangutorum transcriptome and help us to research desert xerophytes’ responses to increased nitrogen in the future.


Introduction
Nitrogen (N) is an essential macronutrient and key signal throughout the plant life cycle [1,2]. In the genetic evolution of plants, the N element is one of the most considerable limiting factors for plant metabolism, as it is a constituent of a variety of many biomolecules, such as amino acids, proteins, chlorophylls, phytohormones, and nucleic acids [3][4][5][6]. N fertilization is known to increase plant yield and productivity [7][8][9][10]; however, over the years, it has also led to increased N fertilizer use by farmers and ranchers. Ecologically, the excessive application of fertilizer has disastrous effects, such as eutrophication [1,11,12], soil acidification [13,14], and air pollution [15,16], as well as changes in the structure and diversity of plant and soil microbials [17,18]. Conversely, N deficiency also affects the activities and processes of plant life by altering the levels of many amino acids and the biosynthesis of some carbohydrates [19,20]. For example, poplar plants reduce N absorption and restrict N metabolism levels under N-deficiency conditions [21,22]. the transcriptome of N. tangutorum, as well as the possible mechanisms that regulate its response to increased levels of N. The information obtained herein may contribute to the exploitation and conservation of various xerophytes.

Growth and Biomass Affected by Nitrogen Treatment
Growth characteristics of N. tangutorum showed differences among the N concentration treatment groups (Figure 1). Compared with the control group (N0), the N6 treatment group had no significant effect on height, basal diameter, leaf biomass, or root biomass of N. tangutorum seedlings. The N36 group showed the highest height, basal diameter, leaf and root biomass, and specific leaf area, their averages being 32.44 ± 1.68 cm, 4.41 ± 0.21 mm, 1.02 ± 0.02 g, 3.30 ± 0.08 g, and 192.18 ± 4.69 cm 2 ·g −1 , respectively. When the N concentration was 60 mmol·L −1 , the root biomass was significantly lower than that of the N0 group, and the biomass was only 56.69% of the N0 group. In addition, the root-shoot ratio of N. tangutorum seedlings decreased gradually with the increase of nitrogen concentration and was significantly lower than that of the N0 group. These results indicated that N36 could be considered as optimal N content for N. tangutorum seedlings growth, N60 as exceeding and depressing.
formation of of N. tangutorum seedlings [53]; in addition, its chlorophyll content and fluorescence activity of leaves were related to the addition of N and P [54]. However, the molecular mechanisms by which N. tangutorum responds to nitrogen addition remains poorly studied. Therefore, in this article, we choose the ideal desert plant, N. tangutorum, as our research material. Our main goal is to share basic information regarding changes in the transcriptome of N. tangutorum, as well as the possible mechanisms that regulate its response to increased levels of N. The information obtained herein may contribute to the exploitation and conservation of various xerophytes.

Growth and Biomass Affected by Nitrogen Treatment
Growth characteristics of N. tangutorum showed differences among the N concentration treatment groups (Figure 1). Compared with the control group (N0), the N6 treatment group had no significant effect on height, basal diameter, leaf biomass, or root biomass of N. tangutorum seedlings. The N36 group showed the highest height, basal diameter, leaf and root biomass, and specific leaf area, their averages being 32.44 ± 1.68 cm, 4.41 ± 0.21 mm, 1.02 ± 0.02 g, 3.30 ± 0.08 g, and 192.18 ± 4.69 cm 2 ·g −1 , respectively. When the N concentration was 60 mmol·L −1 , the root biomass was significantly lower than that of the N0 group, and the biomass was only 56.69% of the N0 group. In addition, the root-shoot ratio of N. tangutorum seedlings decreased gradually with the increase of nitrogen concentration and was significantly lower than that of the N0 group. These results indicated that N36 could be considered as optimal N content for N. tangutorum seedlings growth, N60 as exceeding and depressing. The height and basal diameter growth, (B) leaves and roots biomass, and (C) specific leaf area and root-shoot ratio of N. tangutorum treated with different N additions. Four N addition treatment levels (0 (N0), 6 (N6), 36 (N36), and 60 (N60) mmol·L −1 ) were established in the experiment, with N0 used as the control group. Results are expressed as means ± standard deviation based on three independent experiments. Broken lines with different letters indicate a significant difference (p < 0.05) as determined by analysis of variance and Duncan's multiple range test. The height and basal diameter growth, (B) leaves and roots biomass, and (C) specific leaf area and root-shoot ratio of N. tangutorum treated with different N additions. Four N addition treatment levels (0 (N0), 6 (N6), 36 (N36), and 60 (N60) mmol·L −1 ) were established in the experiment, with N0 used as the control group. Results are expressed as means ± standard deviation based on three independent experiments. Broken lines with different letters indicate a significant difference (p < 0.05) as determined by analysis of variance and Duncan's multiple range test. Table S1 shows the raw reads of N. tangutorum RNA-seq, ranging from 43,516,318 (N36-1) to 57,280,550 (N6-3) under different N addition conditions. The ratio of clean reads to raw reads in each group ranged from 97.10% (N60-2) to 98.49% (N36-3), and the overall average clean reads rate reached 97.93%. These high ratios ensure the splicing of transcripts. Furthermore, the GC content of different treatment ranged from 45.44% (N36-3)
We combined transcript sequences followed by performing hierarchical clustering to obtain 332,420 transcripts and 276,423 unigenes ( Figure 2). We found 86,862 transcripts (31.42% of the total) that were 200-500 bp in length, 77,546 cDNA transcripts that were 500 bp-1 kbp long (28.05% of the total), 71,253 transcripts (25.78% of the total) with a length of 1-2 kbp, and 40,762 transcripts (14.75% of the total) with a length of >2 kbp. In summary, the sequencing and splicing results are of high quality and can be used for the subsequent analysis of gene expression levels. Table S1 shows the raw reads of N. tangutorum RNA-seq, ranging from 43,51 (N36-1) to 57,280,550 (N6-3) under different N addition conditions. The ratio of clean r to raw reads in each group ranged from 97.10% (N60-2) to 98.49% (N36-3), and the ov average clean reads rate reached 97.93%. These high ratios ensure the splicing of scripts. Furthermore, the GC content of different treatment ranged from 45.44% (N to 46.41% (N6-3), with an error rate of 0.03%, which also met our study requirements Supplementary Data Table S1).

RNA-Seq Analysis and Transcript Splicing
We combined transcript sequences followed by performing hierarchical clusteri obtain 332,420 transcripts and 276,423 unigenes ( Figure 2). We found 86,862 transc (31.42% of the total) that were 200-500 bp in length, 77,546 cDNA transcripts that 500 bp-1 kbp long (28.05% of the total), 71,253 transcripts (25.78% of the total) w length of 1-2 kbp, and 40,762 transcripts (14.75% of the total) with a length of >2 kb summary, the sequencing and splicing results are of high quality and can be used fo subsequent analysis of gene expression levels.

GO Functional Annotation
The obtained transcript information was compared with the GO database, an Blast2GO v2.5 software was used to determine the GO classification on all the assem N. tangutorum unigenes. We then identified 122,945 genes that were functionally a tated before subdividing them into 56 GO functional groups. The annotated genes i functional groups were classified according to the three major GO categories (biolo process, cellular component, and molecular function), and the number of funct groups in each major category was 26, 20, and 10, respectively ( Figure 3).

GO Functional Annotation
The obtained transcript information was compared with the GO database, and the Blast2GO v2.5 software was used to determine the GO classification on all the assembled N. tangutorum unigenes. We then identified 122,945 genes that were functionally annotated before subdividing them into 56 GO functional groups. The annotated genes in the functional groups were classified according to the three major GO categories (biological process, cellular component, and molecular function), and the number of functional groups in each major category was 26, 20, and 10, respectively ( Figure 3).

KOG Functional Classification
Comparing the transcript data with information in the KOG database showed that 48,287 unigenes (17.46%) were successfully annotated in the KOG database. According to the KOG functional classification, these genes were arranged into 26 groups (Figure 4). The five groups with the most unigenes, in order of prevalence, were post-translational modification, protein turnover, and chaperones; general functional genes only; translation, ribosome structure, and biogenesis; signal transduction mechanisms; and RNA processing and modification. These contained 6887 (12.8%), 6317 (11.7%), 4801 (8.9%), 3282 (6.1%), and 3273 (6.1%) unigenes, respectively. In addition, 3068 unigenes were related to intracellular transport, secretion, and vesicle control, and 2919 were classified as unknown.

KOG Functional Classification
Comparing the transcript data with information in the KOG database showed that 48,287 unigenes (17.46%) were successfully annotated in the KOG database. According to the KOG functional classification, these genes were arranged into 26 groups ( Figure 4). The five groups with the most unigenes, in order of prevalence, were post-translational modification, protein turnover, and chaperones; general functional genes only; translation, ribosome structure, and biogenesis; signal transduction mechanisms; and RNA processing and modification. These contained 6887 (12.8%), 6317 (11.7%), 4801 (8.9%), 3282 (6.1%), and 3273 (6.1%) unigenes, respectively. In addition, 3068 unigenes were related to intracellular transport, secretion, and vesicle control, and 2919 were classified as unknown.   Figure 5 shows the differences in gene expression among samples treated with different N concentrations. Compared with the N0 group, the number of DEGs in the N6-N36 group grew with the increase of N, and the numbers of upregulated and downregulated genes both expanded with increased N. However, the increase of the number of downregulated genes was significantly larger than that of upregulated genes. In the N60 group, the total number of DEGs was the least; however, the number of upregulated genes was significantly higher than the number of downregulated genes. In addition, principal component analysis (PCA) showed a clear cluster separation of the control (N0) and N treatment (N6, N36, and N60) groups (see Supplementary Data Figure S1).  Figure 5 shows the differences in gene expression among samples treated with different N concentrations. Compared with the N0 group, the number of DEGs in the N6-N36 group grew with the increase of N, and the numbers of upregulated and downregulated genes both expanded with increased N. However, the increase of the number of downregulated genes was significantly larger than that of upregulated genes. In the N60 group, the total number of DEGs was the least; however, the number of upregulated genes was significantly higher than the number of downregulated genes. In addition, principal component analysis (PCA) showed a clear cluster separation of the control (N0) and N treatment (N6, N36, and N60) groups (see Supplementary Data Figure S1). ferent N concentrations. Compared with the N0 group, the number of DEGs in the N6-N36 group grew with the increase of N, and the numbers of upregulated and downregulated genes both expanded with increased N. However, the increase of the number of downregulated genes was significantly larger than that of upregulated genes. In the N60 group, the total number of DEGs was the least; however, the number of upregulated genes was significantly higher than the number of downregulated genes. In addition, principal component analysis (PCA) showed a clear cluster separation of the control (N0) and N treatment (N6, N36, and N60) groups (see Supplementary Data Figure S1).  We used cluster analysis to group similar genes, analyzed the functions of the previously known genes, and predicted those of the previously unknown genes ( Figure 6). Genes clustered in the same group have similar functions and, in some cases, similar expression patterns; they may even participate in the same metabolic processes. As shown in Figure 6, the profiles of gene expression in N6 and N60 were more similar, and they can cluster into the same branch as N0. We used cluster analysis to group similar genes, analyzed the functions of the previously known genes, and predicted those of the previously unknown genes ( Figure 6). Genes clustered in the same group have similar functions and, in some cases, similar expression patterns; they may even participate in the same metabolic processes. As shown in Figure 6, the profiles of gene expression in N6 and N60 were more similar, and they can cluster into the same branch as N0.

Enrichment Analysis of DEGs in KEGG Pathways
The multiple signal transduction and metabolic pathways associated with N. tangutorum DEGs under different N additions were determined through KEGG enrichment analysis. We found that 1101 DEGs (see Supplementary File S1), representing 113 metabolic pathways in the N0 and N6 groups, were annotated in the KEGG database, among

Enrichment Analysis of DEGs in KEGG Pathways
The multiple signal transduction and metabolic pathways associated with N. tangutorum DEGs under different N additions were determined through KEGG enrichment analysis. We found that 1101 DEGs (see Supplementary File S1), representing 113 metabolic pathways in the N0 and N6 groups, were annotated in the KEGG database, among which 14 pathways showed significant enrichment (p-value < 0.05), such as ribosome, glycolysis/gluconeogenesis, and RNA degradation, and so on (Table 1); 2222 DEGs (see Supplementary File S1) were annotated in the KEGG database, representing 121 metabolic pathways in the N0 and N36 groups, and 10 pathways, including porphyrin and chlorophyll metabolism, were significantly enriched. In the comparison between N0 and N60 treatments, 1234 DEGs (see Supplementary File S1) were annotated to 114 metabolic pathways, involving 12 significantly enriched pathways, such as photosynthesis-antenna protein, chlorophyll metabolism, and zeatin synthesis. The significantly enriched KEGG pathways in the other comparison groups are detailed and shown in Table 1.    With protoporphyrin IX as the branch, two metabolic pathways were formed, namely the chlorophyll and haem synthesis pathways. In N6, expression changes in the ferrochelatase   2) were reduced. In N6, the expression of red chlorophyll catabolite reductase (1.3.7.12) and the production of the red chlorophyll-degradation product formed during magnesium removal and transplantbased reactions were reduced, and the synthesis of the primary fluorescence chlorophylldegradation product was inhibited. In addition, the expression levels of chlorophyll synthase (2.5.1.133) and CHIP (1.3.1.111) were downregulated.

Transcriptome Data Verification
We selected 10 DEGs that were highly related to N treatments for qRT-PCR analysis. Figure 8 shows that the RNA-Seq data and expression trends were similar, corroborating the accuracy of the RNA-Seq results, although there were differences in the absolute fold changes between the two methods.

Transcriptome Data Verification
We selected 10 DEGs that were highly related to N treatments for qRT-PCR analysis. Figure 8 shows that the RNA-Seq data and expression trends were similar, corroborating the accuracy of the RNA-Seq results, although there were differences in the absolute fold changes between the two methods.

Discussion
Nitraria tangutorum is the predominant plant class in the Ulan Buh Desert [46,55,56], and its leaves had a high storage capacity for C and N [51]. Previous studies have shown that growth and biomass production of plant was accelerated by N application [57] but was inhibited when the N supply was excessive [24]. Similarly, our study showed that the height, basal diameter, specific leaf area, and leaf biomass of Nitraria tangutorum seedlings were significantly increased by N36, but root biomass and root shoot ratio were inhibited in N60. That is, N36 could be considered as optimal N content for N. tangutorum seedlings growth but more likely to be inhibited under high N concentrations. Moreover, previous studies have shown that nitrogen supplementation not only affects the growth phenotype of plants but also leads to differential expression of their genes. Thus, it is of considerable

Discussion
Nitraria tangutorum is the predominant plant class in the Ulan Buh Desert [46,55,56], and its leaves had a high storage capacity for C and N [51]. Previous studies have shown that growth and biomass production of plant was accelerated by N application [57] but was inhibited when the N supply was excessive [24]. Similarly, our study showed that the height, basal diameter, specific leaf area, and leaf biomass of Nitraria tangutorum seedlings were Plants 2023, 12, 895 11 of 18 significantly increased by N36, but root biomass and root shoot ratio were inhibited in N60. That is, N36 could be considered as optimal N content for N. tangutorum seedlings growth but more likely to be inhibited under high N concentrations. Moreover, previous studies have shown that nitrogen supplementation not only affects the growth phenotype of plants but also leads to differential expression of their genes. Thus, it is of considerable value to explore its gene expression mechanisms in response to increased N under the context of phenotypic differences. In this paper, we studied the transcriptome of N. tangutorum in response to different N additions and established an 89.67 Gb N. tangutorum transcriptome database in which the proportion of clean reads to raw reads reached 97.93%, indicating high sequencing quality [49]. Our study also found that the N36 treatment produced the highest number of clean reads, possibly because N36 induced the diversified expressions of numerous N. tangutorum genes. In addition, more transcripts and unigenes were found in this study of N. tangutorum than in a previous study of Nitraria sibirica under salt stress [58], indicating that N. tangutorum has its own genetic differences in response to different environmental stresses.
In general, plant gene expression is unavoidably affected by a variety of environmental stressors during its growth and development [59,60]. When N availability fluctuates, plants take a number of steps to cope with the new environment [25,61]. For instance, a large number of Populus tomentosa genes were annotated into multiple KEGG pathways, among which multiple pathways related to amino acid and carbohydrate metabolism were significantly enriched, indicating that P. tomentosa shows a significant response to low N stress [20]. In our study, transcriptome data revealed that multiple genes in the leaves of N. tangutorum were differentially expressed under N addition, which was also similar to a previous study that found that the expressions of most ammonia transporter genes in poplar plants were significantly upregulated under low N stress [22]. KEGG metabolicpathway-enrichment analysis for DEGs showed that the porphyrin/chlorophyll metabolic pathways were significantly enriched following each treatment; thus, we speculated that these responses may be associated with the diversity of enzymes (e.g., GS) in plants, reflecting the complexity of their roles in plants growth [24,62]. Other studies have shown that the plant "ribosome" pathway has undergone considerable changes under abiotic stress (e.g., drought stress [63]). A similar phenomenon was reflected in this study; the ribosome pathway was the KEGG pathway with the most significant difference between N0 and N6, indicating that a small N addition was more conductive to the differential expression of N. tangutorum ribosomal-pathway-related genes. Moreover, Anthocyanins are natural colorants belonging to the flavonoid family that have been shown to possess potent antioxidant properties [64]. Wang et al. [52] found that the "lavonoid synthesis" pathway was significantly enriched, and an anthocyanin synthesis gene, Oxoglutarate/iron-dependent dioxygenase (2-GO), was annotated by GO and highly upregulated, and similar to our study, they also found the anthocyanin synthesis, flavonoid biosynthesis, photosyntheticantenna protein, amino acid biosynthesis, and metabolism processes were also significantly enriched. Thus, we speculate that upregulation of anthocyanin-related genes indicated that anthocyanins play an important role in reactive oxide species scavenging in N. tangutorum under N addition. In a word, the addition of N leads to changes in the external growth environment of N. tangutorum, which may make N. tangutorum more sensitive to external environmental stress and ultimately lead to a more frequent adjustment of transcriptional output [65].
Chlorophyll is the key factor involved in plant photosynthesis, providing energy for plant growth, development, and productivity [66][67][68]. As a component of the chloroplasts, increased nitrogen is beneficial for chlorophyll synthesis, up to a point [69][70][71]. As is well known, the molecular regulation of chlorophyll biosynthesis is a complex process, affected by the external abiotic factor and regulated by related genes [72,73]. The first step in chlorophyll biosynthesis is ALA synthesis, which is synthesized from glutamic acid tRNA synthetase, δ-ketoglutaric acid tRNA synthetase, GluTR tRNA reductase, and GSA-AT. In this paper, our data analysis showed that compared with the control, GluTR was not expressed in N6 but induced in N36 and was further upregulated in N60. With the increase in N addition, the expression trend of glutamyl tRNA reductase was consistent with glutamine-1-hemialdehyde transaminase. These results were not completely consistent with the expression of genes related to chlorophyll biosynthesis in maize leaves under zinc stress [74], indicating that different abiotic stresses had different effects on chlorophyll synthesis in plants. Thus, we predicted that the expression levels of these two enzyme genes in N. tangutorum would increase with increased N addition. Previous reports have shown that the expression of ALA-synthesis-related genes can affect the chlorophyll contents of plants [41] and that GluTR can regulate ALA synthesis at the transcriptional level, thereby affecting chlorophyll synthesis [75,76]. Here, we obtained similar results and speculated that glutamyl tRNA reductase is likely to be a key regulatory site for ALA synthesis when N. tangutorum receive external N input.
The second considerable pathway for chlorophyll synthesis involves the synthesis of protoporphyrin IX [77]. This pathway begins with ALA, from which protoporphyrin IX is formed through isomerization, decarboxylation, and oxidation reactions catalyzed by bile pigment synthase, bile pigment deaminase, uroporphyrin proporphyrin III synthase, uroporphyrin proporphyrin III decarboxylase, coproporphyrin III oxidase, and protoporphyrinogen III oxidase. In the present study, the bile pigment synthase gene was not expressed in N6 but was low-expressed in N36 and overexpressed in N60, indicating that the N concentration in the environment determined the bile pigment synthase expression level in N. tangutorum leaves. In addition, the gene-encoding bile chromatogen deaminase was downregulated in N6 and N36 but not expressed in N60. Both up/downregulation of the uroporphyrin proporphyrin III decarboxylase gene were observed in N36; however, there was more downregulation than upregulation. Therefore, we speculate that upregulation of the bile pigment deaminase and uroporphyrin proporphyrin III decarboxylase genes could occur at appropriate levels following N addition. However, the specific regulation mode of these differentially expressed genes remains to be further studied.
The coproporphyrin III oxidase gene was downregulated in N6 and N36, and ProIX oxidase was upregulated in N60; thus, it is possible that these two genes coordinate with each other before ProIX synthesis to increase ProIX expression under high N concentration. When ProIX production is catalyzed by ferrochelatase, it forms ferrous heme and then enters the heme synthesis pathway; when ProIX encounters magnesium chelase, Mg-protoporphyrin IX enters the chlorophyll synthesis pathway [78]. Previous findings showed that plants controlled the flow of ProIX towards chlorophyll synthesis by regulating magnesium chelase [79]. Similarly, the affinity of iron chelase for ProIX is lower than that of magnesium chelase, indicating that most ProIX leads to chlorophyll synthesis [80]. Furthermore, these genes were involved mainly in the chlorophyll synthesis pathway in the three experimental groups, and the expression of genes related to heme synthesis was very low. Compared with the controlled plants, the magnesium chelatase H subgroup gene of the N6 group was downregulated. Under increasing N, the upregulated expression of this gene was observed in N36, and although both up/downregulation of this gene were observed in N60, it was less upregulated than in N36, which led us to speculate that the expression of a magnesium-chelating enzyme during N processing is beneficial for chlorophyll synthesis over heme synthesis. Similarly, Liu et al. [41] found in the study of Brassica that the CAB domain in BrFC2 was not the structure to maintain the catalytic activity of heme enzyme, and only after the single base mutation of dBrFC2 could both increase the content of chlorophyll and heme in plants. It can be seen that N36 treatment is likely to change the domain of one or more genes in the chlorophyll synthesis of N. tangutorum but did not cause mutations in genes related to heme synthesis. Previously, it was shown that chlorophyllin a synthesis was catalyzed by magnesium chelase and other enzymes, such as POR [81]. Here, we found that the expression of protochlorophyllide reductase was lower in both N6 and N60 and that the degree of downregulation was almost the same in each case, whereas the enzyme was not expressed in N36. These results indicated that original chlorophyll ester reductase expression affected chlorophyll ester a synthesis by reducing its expression at both low and high nitrogen concentrations. In addition, the translation of chlorophyll b into chlorophyll a is likely to be part of the chlorophyll degradation pathway [82]. In the present study, relative to the control, we found that chlorophyll b reductase gene expression was lower in N6, N36, and N60, and the capability to translate chlorophyll b to chlorophyll a was reduced at higher N concentrations. Thus, we speculated that increased nitrogen addition could reduce chlorophyll degradation.  [47]. To reduce the impact of the experimental seedling growth difference on the research results, uniformly growing seedlings were selected from the above-bred half-sib families, then transplanted into self-developed PVC barrels (40 cm height × 16 cm in diameter) on 8 May 2015, with one plant per barrel. The soil matrix (the main physicochemical properties of the soil matrix are given in Table 2) consisted of local abandoned cornfield topsoil and river, which was sand mixed 1:1 (v/v) and screened. On 23 May 2015 (15 days after transplanting), the seedlings in the PVC barrels were treated with N fertilizer (N content 46.7% of urea (Shaanhua Coal Chemical Industry Group Co., Ltd., Shanxi, China)). Four additional N treatment levels (0 (N0), 6 (N6), 36 (N36), and 60 (N60) mmol·L −1 ) were established in the experiment, each with 25 replicates per treatment, with N0 used as the control group. To prevent the effects of natural precipitation on the experimental results, a transparent canopy was built and used to control the water supply. On 1 August 2015 (70 days after N treatment), leaf samples were collected from each treated seedling between 9:30 and 10:30 am. We randomly selected 18 barrels of seedlings from each treatment, 6 mature leaves were collected from each seedling, and a total of 108 leaves were collected for each treatment. We divided the 6 leaves of each seedling into 3 two-leaf repetitions. Each treatment therefore had three repetitions with 36 mixed leaves per repetition. Then, the samples were quickly stored at −80 • C until use. In addition, we measured and counted the plant height, basal diameter, specific leaf area, leaf biomass, root biomass and root-shoot ratio of N. tangutorum seedlings with different N treatments before sampling.

Library Preparation for Transcriptome Sequencing
Total RNA of N. tangutorum was isolated from their leaf powder (three biological plants replicated per treatment) using the RNAsimple Total RNA Kit (Tiangen, Beijing, China). Then, 1.5 µg RNA was taken from each sample as the input material for the RNA preparations. Sequencing libraries were generated according to Liu et al.'s method [47]. The mRNA was isolated and purified from the total RNA using poly-toligo-attached magnetic beads. Next, first-strand cDNA was synthesized using reverse transcriptase, and second-strand cDNA was synthesized using DNA polymerase I and RNase H (Novogene Technology Co., Ltd., Beijing, China). The library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, MA, USA) to obtain the target cDNA fragments. Then, 3.0 µL USER Enzyme (NEB, Ipswich, MA, USA) was used with sizeselected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR.
Next, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primer, and Index (X) Primer. Library quality was evaluated on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) and, finally, PCR amplification was performed to obtain the complete cDNA library. All of these experiments were performed with the help of Novogene Technology Co., Ltd. (Beijing, China).

Enrichment Analysis of Differentially Expressed Genes
The analysis and screening of the differentially expressed genes (DEGs) were performed using the DESeq software [83]. The numbers of DEGs between the N0 groups and N addition treatment groups were counted, including the up/downregulated genes. The p-value was adjusted using the Benjamini and Hochberg method. The corrected p-value of <0.01 and a log2-fold change > 1 were set as the threshold for significant differential expression [47]. ImageGP (http://www.ehbio.com/ImageGP, accessed on 16 October 2015) was used to achieve the GO analyses [29]. All the obtained DEG sequences were annotated in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database by KOBAS (2.0) (http://kobas.cbi.pku.edu.cn/, accessed on 23 November 2015) to identify the signal transduction and metabolic pathways involved in DEGs. When the p-value threshold of ≤0.05 was obtained, the KEGG analysis was considered significantly enriched by DEGs [84]. All raw data obtained from the experiments described above are available for review in Supplementary File S2.

Verification of RNA Sequencing Data by RT-qPCR
To verify the RNA-Sequencing data, 10 candidate DEGs (1, 7, and 2 DEGs from the N0 vs. N6, N0 vs. N36, and N0 vs. N60 comparison groups, respectively) of N. tangutorum were randomly selected, and qRT-PCR analysis was conducted to validate the differences in their expression levels. The specific primers were designed using PRIMER5 software (http://www.PremierBiosoft.com, accessed on 17 January 2016) [85] and are listed in Supplementary Data Table S2. We reverse-transcribed 1.0 µg of the total RNA of each N. tangutorum sample with the Goldenstar TM RT6 cDNA Synthesis Kit (Novogene Technology Co., Ltd., Beijing, China); then, the cDNA was amplificated using 2 × T5 Fast qPCR Mix (SYBR Green I). PCR amplification was performed under the following conditions: one cycle of 1 min at 95 • C, followed by 40 cycles at 95 • C for 10 s, 60 • C for 5 s, and 72 • C for 10 s. Afterward, the relative gene expression levels were analyzed using the 2 −∆∆Ct method [86]. The results were analyzed by using the Bio-Rad CFX Manager 3.1 software (Bio-Rad, Hercules, CA, USA). qRT-PCR for each N. tangutorum sample was repeated three times.

Conclusions
A total of 14,170 DEGs were identified. There were 1101, 2222, and 1234 DEGs assigned to the N0 vs. N6, N0 vs. N36, and N0 vs. N60 groups in the 113, 121, and 114 metabolic pathways classified in the KEGG database, respectively. The metabolic pathways upregulated by increased nitrogen included anthocyanin biosynthesis, carotenoid biosynthesis, porphyrin and chlorophyll metabolism, flavonoid biosynthesis, and amino acid metabolism. The magnesium chelatase H subunit was upregulated in N36, promoting chlorophyll a synthesis, implying that N36 metabolism favored chlorophyll synthesis over heme synthesis.